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ABSTRACT 

We analyze the formation and evolution of stellar bars in galactic disks embedded in mildly triaxial 
cold dark matter (CDM) halos that have density distributions ranging from large flat cores to cuspy 
profiles. We have applied tailored numerical simulations of analytical and live halos which include the 
feedback from disk/bar system onto the halo in order to test and extend earlier work by El-Zant & 
Shlosman (2002). The latter employed the method of Liapunov exponents to analyze the fate of bars 
in analytical asymmetric halos. We find the following: (1) The bar growth is very similar in all rigid 
axisymmctric and triaxial halos. (2) Bars in live models experience vertical buckling instability and 
the formation of a pseudo-bulge with a boxy/peanut shape, while bars in rigid halos do not buckle. 
(3) In live axisymmetric halos, the bar strength varies by a factor of <j 2, in growth or decay, during 
the secular evolution following the buckling. The bar pattern speed evolution (i.e., deceleration) 
anticorrelates with the halo core size. In such halos, the bar strength is larger for smaller disk-to-halo 
mass ratios D/H within disk radii, the bar size correlates with the halo core sizes, and the bar pattern 
speeds correlate with the halo central mass concentration. In contrast, bars embedded in live triaxial 
halos have a starkly different fate: they dissolve on a timescale of ^ 1.5 — 5 Gyr due to the onset of 
chaos over continuous zones, sometimes leaving behind a weak oval distortion. The onset of chaos is 
related to the halo triaxiality, the fast rotating bar and the halo cuspiness. Before the bar dissolves, 
the region outside it develops strong spiral structures, especially in the live triaxial halos. (4) More 
angular momentum is absorbed (fractionally) by the triaxial halos as compared to the axisymmetric 
models. The disk-halo angular momentum exchange is mediated by the lower resonances in the 
latter models. (5) Cuspy halos are more susceptible than flat-core halos to having their prolateness 
washed out by the action of the bar. The subsequent evolution is then similar to the case of a cuspy 
axisymmetric halos. We analyze the above results on disk and bar evolution in terms of the stability of 
trajectories and development of chaos in the system. We set important constraints on the triaxiality 
of DM halos by comparing our predictions to recent observational results on the properties of bars 
out to intermediate redshifts z ~ 1. 

Subject headings: galaxies: evolution - galaxies: formation - galaxies: halos - galaxies: kinematics 
and dynamics - galaxies: structure - cosmology: dark matter 



1. INTRODUCTION 

Stellar bars are recognized as the single most impor- 
tant internal factor which drives the evolution of disk 
galaxies both dynamically and secularly, modifying their 
morphology in this process. Modern understanding of 
the bar growth is based on the efficiency of angular mo- 
mentum exchange between the (bar-forming) inner disk 
and the surrounding dark matter halo, outer disk, bulge, 
and, to a certain degree, with the immediate galactic 
environment (Athanassoula 2003). A number of yet to 
be identified and investigated intricacies can affect the 
efficiency of this process. 

In this work we attempt to analyze the formation and 
evolution of stellar bars embedded in fully grown disks 
in the presence of triaxial ^ halos with various radial den- 
sity profiles, from cuspy to flat-core models. We study 
both analytical (i.e., rigid) and live triaxial halos and 
compare our results with the evolution of the bars in ax- 

^ We also use the term 'asymmetric' concurrently with 'triaxial' 



isymmetric ones. Corollaries for disk evolution, such as 
the back-reaction of disks and bars on the halo shapes, 
are investigated as well. This effect has broad implica- 
tions for the cosmological evolution of galaxies. Addi- 
tional evolutionary effects on the 3-D structure of stellar 
bars embedded in such halos will be addressed elsewhere. 

The evolution of bars is expected to be substan- 
tially altered if the surrounding halos are even mildly 
non-axisymmetric (El-Zant & Shlosman 2002, hereafter 
ES02). However, this effect in a live disk- halo system was 
not verified so far. Competing gravitational torques from 
a bar and a halo acting on the main families of planar 
and 3-D periodic orbits can destabilize them and dra- 
matically reduce their ability to trap neighboring trajec- 
tories, thus inducing chaos and dissolving otherwise sta- 
ble structural features in disk galaxies. The full extent 
of these non-linear effects is not yet clear, but one can 
expect them to speed up secular changes in the collision- 
less components and to facilitate the angular momentum 
loss in the gaseous one. Although triaxial shapes of dark 
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matter halos appear inherent in the cosmological numer- 
ical simulations (e.g., Barnes & Efstathiou 1987; Frank 
et al. 1988; Dubinski & Carlberg 1991; Warren et al. 
1992; Cole & Lacey 1996), relatively little attention has 
been paid so far to this issue when building self-consistent 
galactic models. 

Observational constraints on the shapes of galactic ha- 
los coming from gas kinematics in disk galaxies (e.g., 
Sparke 1986; Sackett 1999; Andersen et al. 2001; Merri- 
field 2002; Dckel & Shlosman 1983) and their polar rings 
(Sackett et al. 1994; Sackett & Pogge 1995), gravita- 
tional lensing (Kochanek 1995; Oguri et al. 2003; Hoek- 
stra, Yee & Gladders 2004) and X-ray gas in ellipticals 
(Buote & Canizares 1994, 1996; Buote et al. 2002) re- 
main inconclusive. 

Nevertheless, some clues exist. Residual potential ax- 
ial ratios (both flatness and prolateness^) of about 0.9 in 
Cold Dark Matter (hereafter CDM) halos are plausible, 
even in present day galaxies (e.g., Kuijken & Tremaine 
1994 [for the Milky Way]; Rix & Zaritsky 1995; Helmi 
2004). Another evidence for triaxiality of the DM halo 
is provided by the X-ray isophote position angle twist 
observed by Chandra in NGC 720 (but not in stellar 
isophotes!) (Buote et al. 2002). A nearly prolate DM 
halo, with T ~ 0.985 (and b/a - 0.791, c/a ~ 0.787) has 
been inferred for the elliptical galaxy NGC 5128 from a 
kinematic study of planetary nebula system (Peng et al. 
2004). We conjecture, therefore, that while some individ- 
ual halo shapes are compatible with being mildly prolate 
even at the present time, the statistical significance of 
this effect is not clear, and it is not clear whether disk 
galaxies differ in this respect from ellipticals. Overall, 
the prolateness of contemporary halos appears to be in- 
significant. 

The main motivation behind this work is that theo- 
retically a galactic halo is expected to acquire its triax- 
ial shape during its initial collapse and to support this 
shape during the ongoing process of a hierarchical merg- 
ing. The degree of triaxiality'^ will depend on the merging 
history, specifically on the relative angular momenta of 
the merger precursors and interaction frequency — which 
is difficult to quantify. 

Numerical simulations indicate that at very high red- 
shifts, at the epoch of galaxy formation, the halos can 
be significantly triaxial. Detailed properties of individ- 
ual numerical halos, such as triaxial shapes and radial 
density structure, can be extracted from the models and 
directly confronted by their observational counterparts. 
Flatness and prolateness in the models appear to increase 
with the halo mass, being somewhat milder in the ACDM 
than CDM cosmology (e.g., for a recent review Bullock 
2002). This triaxiality seems to be marginally higher in 
the outer halo parts, independently of the halo mass and 
of the ratio of rotational-to-dispersion velocity, indicat- 
ing that the halos are supported by anisotropic velocity 
dispersions. 

The high halo triaxiality that may be present at the 
epoch of early galaxy formation can be substantially di- 

^ We define the halo flatness as / = 1 — c/a and its prolateness 
(i.e., equatorial ellipticity) as en = I — b/a 

3 The triaxiality is defined here as T = [1 - (fe/a)2]/[l - (c/a)^], 
where c/a is halo's polar-to-longest equatorial axis ratio and b/a - 
equatorial axis ratio. T = 1 corresponds to a prolate halo, while 
T = to an oblate one 



luted by the present-day. For instance, the addition of 
spherical and/or axisymmetric baryonic components to 
the system (e.g., during the formation and development 
of a galactic disk) can wash out the halo triaxiality, still 
keeping it non-negligible (Dubinski 1994; Kazantzidis et 
al. 2004). 

From an observational standpoint it appears that pro- 
files of the DM halos in fully formed galaxies tend to 
have nearly constant density cores (Flores & Primack 
1994; Moore 1994; Burkert 1995; Kravtsov et al. 1998; 
Borriello & Salucci 2001; de Blok & Bosma 2002; Gen- 
tile et al. 2004; etc.). A sim ilar effect has been recently 
observed in galactic clusters (Sand et al. 2002). At the 
same time, dissipationless CDM simulations of galactic 
halos agree with a universal density profile oc r~", where 
a = 1 — 1.5 (e.g., Dubinski & Carlberg 1991; Warren et 
al. 1992; Crone et al. 1994; Cole & Lacey 1996; Tormen 
et al. 1997; Huss et al. 1999; Fukushige & Makino 1997; 
Moore et al. 1999; Klypin et al. 2000; Jing & Suto 2002; 
Power et al. 2003). (Note, Jing & Suto (2000) claim that 
the profiles are not universal.) Navarro, Frenk & White 
(1996, hereafter NFW) have found a fitting formula for 
the density profile of DM halos, for a wide range of cos- 
mological models, in which the inner profile diverges as 
r~^, while the outer profile drops as r~^. It has been 
also demonstrated theoretically that a cuspy density pro- 
file arises inherently from the cold gravitational collapse 
in an expanding universe (Lokas & Hoffman 2000). The 
CDM model, therefore, predicts that the inner density 
profile of galactic scale DM halos is characterized by a 
density cusp while observations of the dynamics of the 
central regions of galaxies imply a core-halo structure of 
the DM. Another disagreement with observations is the 
so-called angular momentum catastrophe — the iV-body 
and gas dynamical simulations consistently result in too 
small galactic disks due to the overall lack of angular mo- 
mentum necessary to reproduce the observed disk sizes 
(e.g., Burkert & D'Onglia 2004). 

This controversy between observations of DM cores 
and density cusps in numerical models is not a funda- 
mental one and can be resolved within the general con- 
text of CDM cosmology. Within the conventional physics 
framework, interactions with the dissipative and clumpy 
baryonic component, such as dynamical friction, dur- 
ing the initial stages of collapse can level off the central 
density cusps (which have been shown to be therniody- 
namically improbable) and produce harmonic cores in 
DM halos (El-Zant, Shlosman & Hoffman 2001). Even 
though this may affect the isodensity contours, making 
them rounder, it need not symmetrize the isopotentials; 
these can remain asymmetric if triaxiality is not affected 
beyond some radius (for example, a uniform bar has a 
non-axisymmetric force contribution inside its density 
figure). More generally, this process was shown to re- 
place the DM cusps with baryonic cusps (El-Zant et al. 
2004). Furthermore, asymmetric and flat core halos can 
have interesting implications for the disk growth and cor- 
relate the properties of the central supermassive black 
holes with those of galactic bulges and halos themselves 
(El-Zant et al. 2003). Alternatively, the central density 
cusps have been proposed to dissolve by the action of 
galactic bars (Weinberg & Katz 2002; but see McMillan 
& Dehnen 2005). 

A number of approaches can be taken in order to 
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construct triaxial halos and elliptical galaxies (e.g., 
Aarseth & Binney 1974; Merritt & Valluri 1996; HoUey- 
Bockelmann et al. 2001; Moore et al. 2004) and to study 
the bar dynamics within triaxial halos. Ideta & Hozumi 
(2000) have used analytical density distributions for two 
highly prolate halos with (equatorial) axial ratios of 0.6 
and 0.75 and a very steep density distribution outside 
the core. Here we choose to construct stable triaxial 
models with a subsequent introduction of axisymmetric 
collisionless disks which are then released and their time 
evolution is observed. In our initial conditions we tend to 
follow ES02, who studied systematically the stability of 
stellar bars in analytical halos of different triaxiality and 
central concentration by means of Liapunov exponents. 
While in ES02 approach the feedback from disk/bar sys- 
tem onto the halo is naturally ignored, our models pre- 
sented here fully account for it. 

In ES02, mildly triaxial shapes have been used, with 
the gravitational potential axial ratios in the equatorial 
plane, b/a, between 0.9 and 1.0, and c/a = 0.8, where c is 
the halo polar axis. A clear trend has been found, in the 
sense of models becoming intrinsically chaotic with grow- 
ing triaxiality and central concentration. For small halo 
(flat) core sizes, ~ 0.5 — 2 kpc, and potential axis ratios 
of 0.9 — 0.95 most of the trajectories integrated appeared 
chaotic and had large Liapunov exponents. Importantly, 
trapping by neighboring stability islands is insignificant, 
because the distribution of values of Lyapunov exponents 
is very similar to the distribution of occupied configura- 
tion space volume. This means that stellar bars under 
these conditions would disintegrate on timescales of a few 
dynamical times, much shorter than the Hubble time, 
as chaotic trajectories quickly diffuse out of the bar's 
configuration space. Even spherically-symmetrical mod- 
els, with a small core size, showed a healthy fraction of 
chaotic trajectories, though connected regions of regular 
orbits aligned with the bar remained in this case and a 
self-consistent bar could be maintained. 

What are the reasons for a dramatic increase in the 
fraction of chaotic orbits in the barred disks with the 
increase in central concentration (i.e., cuspiness) and tri- 
axiality in the halo models of ES02? First, in terms of in- 
variants of motion, a stable 3-D figure must be built from 
trajectories which at least approximately conserve them. 
The chaos appears when the number of invariants of mo- 
tion becomes smaller that the dimensionality of the sys- 
tem. While in the flat core systems, the potential can be 
approximated as quadratic (i.e., harmonic) and motions 
along the coordinates are independent of each other, in 
cuspy potentials these motions are coupled, which leads 
to the overall decrease in the number of invariants of mo- 
tion and typically to a chaotic behavior. In other words, 
cuspy density distribution cause solutions for the Pois- 
son equation to be far from quadratic, i.e., when the 
potential is expanded in power series, it will have a non- 
negligible contribution from terms beyond the quadratic 
one. These terms produce the coupling between different 
degrees of freedom in equations of motion which become 
highly non-linear. Such systems are prime candidates to 
excite chaotic motions when perturbed. 

The second reason is related to the time-dependent 
character of the azimuthal force field comprised from the 



fast rotating^ bar and a non- or slowly tumbling halo. In 
other words, the origin of chaos in this case lies in the 
comparable (in value) and competing forces from the bar 
and the asymmetric halo. (Note, that in a halo potential 
only, most of the orbits are regular, even in the triaxial 
case.) 

This onset of chaos in the presence of a bar within 
even a mildly triaxial halo hints that such configura- 
tions are structurally unstable on dynamical timescales. 
Since numerical cosmological halos arc both triaxial and 
centrally-concentrated, serious questions arise about the 
survival of large-scale stellar bars or of the halos 's tri- 
axiality under these conditions. This issue defines the 
general thrust of our work. 

This paper is structured as follows. In the next section 
we provide the initial conditions for numerical modeling. 
Sections 3 and 4 describe our results for analytical and 
live models, and section 5 is devoted to discussion and 
concluding remarks. 

2. NUMERICAL MODELING 

We have introduced the following dimensionless model 
units. The spatial distance unit is taken as r = 10 kpc, the 
mass unit is M = 10^^ Mq and the gravitational constant 
is chosen to be G = 1, which result in a time unit of 
T= (r'^/GM)^/^ = 4.7 X 10''yrs, corresponding to the 
dynamical timescale, r^yn- In these units, the velocity is 
given in 208 km s~^. The actual physical units are used 
when it is needed for clarity. 

We have used version FTM-4.4 of A^-body code (e.g.. 
Heller & Shlosman 1994; Heller 1995) with N = 6-9x10^ 
particles and gravitational softening of 100 pc to simu- 
late the collisionless disk and spheroidal galactic compo- 
nents (i.e., bulges and DM halos) in a large number of 
models. Our results appear to be reasonably indepen- 
dent of A^. The gravitational forces have been computed 
using Dehnen's (2002) falcON force solver, a tree code 
with mutual cell-cell interactions and complexity 0(N). 
It conserves momentum exactly and is about 10 times 
faster than optimally coded Barnes & Hut (1986) tree 
code. 

To analyze the angular momentum redistribution in 
the models, we have applied the spectral analysis method 
described in Binney & Spergel (1982), in conjunction 
with our nonlinear orbit finding algorithm (e.g., Heller 
& Shlosman 1996). Athanassoula (2002) has shown that 
the angular momentum exchange between the disk and 
the halo is mediated by the lower resonances, with the 
disk losing and the halo gaining the angular momentum. 
Overall, the resonances provide an independent test of 
the quality of the numerical scheme. When the A^-body 
potential is too 'grainy', the particles cannot be locked 
by the resonances and so the latter efficiency in redis- 
tributing the angular momentum is sharply reduced. 

2.1. Initial conditions 

The main challenge in running evolutionary models of 
galactic disks embedded in triaxial halos is to form the 
initially stable halo configurations in the first place. A 
limited number of options known to us was not consid- 
ered to be satisfactory, as they did not allow for a suf- 

By 'fast rotating bar' we mean a bar which extends to about 
its corotation radius 




Fig. 1. — Top panels: Density distribution of the live triaxial halos along its principal axes (solid line: along a, dashed line: along 6, 
dotted line: along c-axis) for models with a varying core size, Rn. Bottom panels: Potential axes ratios /3 = b/a (solid line) and 7 = c/a 
(dotted line) for the same models. 




Fig. 2. — Initial rotation curves for {upper panel:) RSI, RS4 and RS5 models (rigid spherically-symmetric halos, constant A/d/A/h) 
within r = 10 kpc; {middle panel:) RSI to RS3 models (rigid spherically-symmetric halos, constant Vh); and {lower panel:) LSI to LS 3 
(live axisymmetric halos) with cores of 10, 2 and 0.5 kpc. Thin solid line shows halo contribution, dotted line — disk contribution, and 
thick solid line — the total. The units of length and velocity correspond to our model units. 



ficient control of initial triaxiality and mass distribution 
in the halos. Instead we have designed the following pro- 
cedure to obtain the required halo properties, which is 
described below. 

First, we have used the known density-potential pairs 
of required symmetry to lay out particles up to some 
given radial cut-off radius. The velocity distribution 
function has been taken oc E~°', where E is the energy 
and a = 2.5 — 3.5, i.e., below the 3.5 value for the Plum- 



mer sphere. Next, the halo has been evolved for about 
50 Tdyn to allow for any residual relaxation to cease. This 
has resulted in a stable spherically-symmetric configura- 
tion. For the axisymmetric halos, a frozen stellar disk 
has been gradually introduced thereafter and the halo 
was been given time to relax again in this growing disk 
potential. The disk potential slightly changes the initial 
halo core size and therefore we apply an adiabatic cool- 
ing procedure to the halo particles within the core region. 
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TABLE 1 
INITIAL MODEL PARAMETERS 



iVioacl 


Hu [kpc\ 




P 


7 


IN OtCS 


RS 1 


10.0 


1.0 


1.0 


1.0 


Rieid sphericallv-svmmetric halo 


RS 2 


2.0 


1.0 


1.0 


1.0 


Rigid spherically-symmetric halo 


RS3 


0.5 


1.0 


1.0 


1.0 


Rigid sphericallv-svmmetric halo 


RS 4 


2.0 


0.721 


1.0 


1.0 


Rigid spherically-symmetric halo 


RS 5 


0.5 


0.708 


1.0 


1.0 


Rigid spherically-symmetric halo 


LS 1 


10.1 


1.13 


1.0 


0.99 


Live axisymmetric halo 


LS 2 


1.8 


1.18 


1.0 


0.98 


Live axisymmetric halo 


LS3 


0.3 


1.20 


1.0 


0.98 


Live axisymmetric halo 


RTl 


10.0 


1.0 


0.9 


0.9 


Rigid triaxial halo 


RT2 


2.0 


1.0 


0.9 


0.9 


Rigid triaxial halo 


RT3 


0.5 


1.0 


0.9 


0.9 


Rigid triaxial halo 


RT4 


2.0 


0.719 


0.9 


0.9 


Rigid triaxial halo 


RT5 


0.5 


0.708 


0.9 


0.9 


Rigid triaxial halo 


LTl 


10.5 


1.14 


0.79 


0.72 


Live triaxial halo 


LT2 


2.4 


1.19 


0.84 


0.77 


Live triaxial halo 


LT3 


0.5 


1.19 


0.88 


0.79 


Live triaxial halo 


LT3HM 


0.7 


1.10 


0.84 


0.77 


as LT 3, but with half-mass disk 


LT3MT 


0.5 


0.4 


0.9 


0.8 


as LT 3, but maintaining triaxiality 



Note. — Columns: (1) model type (see text); (2) (fitted) halo core radius; (3) (fitted) constant in eq. 1; (4) (fitted) halo equatorial 
axial ratios; (5) (fitted) halo polar-to-longest axis ratios 



by gradually reducing the velocities. Before introducing 
the disk for the triaxial cases, we have implemented a 
similar adiabatic 'heating-cooling' procedure to the halo 
particles to create a triaxial halo of required flatness and 
prolateness, by using 'heating' of the velocities along the 
X-axis and 'cooling' along both the y- and the z-axes, 
over some period of time, i.e. some 40 Tdyn- The frac- 
tional energy heating/cooling per dynamical time is very 
small, ^ 10^^ per machine timestep (corresponding to 
about 10~^ per dynamical time), so the system is never 
taken out of a dynamical equilibrium. The disk poten- 
tial partially dilutes the halo triaxiality and therefore wc 
apply a second heating-cooling procedure after the disc 
has been introduced in the models. This finally results 
in both the required (3 = b/a and 7 = c/a distributions 
and core sizes in the halo potential (Fig. Note that 
in the models with smaller cores it is more difficult to 
impose (and maintain) the triaxiality at all radii. 

For a direct comparison with ES02 in this work, we 
aimed to closely reproduce their initial conditions. A 3- 
D mass distribution (live or analytical) which pairs with 
the nonrotating logarithmic potential, <I>h, (e.g., Binney 
& Tremaine 1987) was used for the halo, 

*H = ^ log {Rl + x^ + p-^v^ + j-^z^) , (1) 

where Vjj is the asymptotic (in the limit R ^ i?H and 
/3,7 1) circular velocity, i?H is the halo core radius 
and P (= b/a), 7 (= c/a) < 1 are the azimuthal and 
polar potential axis ratios. Note that the correspond- 
ing mass distribution is diverging for this potential. We, 
therefore, applied a truncation radius for models with 
live halos, i.e. 50 and 100 kpc for the axisymmetric and 
triaxial halos, respectively. Although the potential of 
these models is going to change, both as a result of the 
disk addition and, in the case of live halo response, to 
the disk evolution, we find it nevertheless beneficial to 
compare the evolving potential to the fitted logarithmic 
one. For spherical (axisymmetric) models P ~ j ~ 1.0, 



while for triaxial models P = 0.9 and 7 = 0.8, approxi- 
mately for live models, after the halo has relaxed in the 
disk potential. 

We have used two sets of initial halo models. In the 
first set, 1 ^ 4 ^ 5 (see Table 1) — for rigid halos only 
- we adjust Vh requiring 1 : 1 ratio of halo-to-disk mass 
within r ~ 1. This means that on the scale of r 1 
models with different core size have the same mass con- 
centration. On smaller scales of course the more cuspy 
models will be always more centrally concentrated, which 
can affect the properties of developing stellar bars. 

In the second set, 1 ^ 2 — > 3, for which we run both 
rigid and live models, we retain the value of Vn, obtained 
for an axisymmetric halo with a core radius of i?H = 
1. For further halo configurations, the total halo mass 
within its truncation radius r = 100 kpc, is 10.0 in our 
units. Thus, when moving from a large flat core to cuspy 
models, we move from a maximal disk model to a halo- 
dominated one. 

Similarly, as in ES02, the stellar disk has been set up 
following the potential in the form of (Miyamoto & Nagai 
1975), 



which describes a disk-bulge system, with the parameters 
Ad = 0.284 and Bd = 0.05, determining the scalelength 
and height, respectively. The disk mass is 0.5 within 
r = 1 and about 0.6 within the disk cut-off radius, i.e. r = 
2.5. The initial Toomre's Q parameter is taken 1.2 and 
constant with radius. The disk rotational velocities in 
triaxial halos have been set using our standard procedure 
but ignoring the mild triaxialities in the total potential. 
The initial models are summarized in Table ^ a-nd the 
initial rotation curves are shown in Fig. |2| 
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3. RESULTS 

We now present the bar evolution and its interaction 
with the outer disk and halo for the numerical models 
listed in Tabled We start with disks in rigid axisymmct- 
ric^ halos and compare them with those in rigid triaxial 
halos, assuming constant disk-to-halo mass ratios within 
10 kpc, and varying the flat core sizes from large 10 kpc 
to cuspy ones. We then continue with live axisymmetric 
and triaxial halos. 

3.1. Comparing Bars in Rigid Axisymmetric and 
Triaxial Halos with Constant Disk-to- Halo Mass 
Ratios 

3.1.1. Rigid axisymmetric halos 

Only circumstantial evidence exists comparing the 
bar formation and evolution in analytical (i.e., rigid) 
and live halos, even in the case of an axial symmetry 
(e.g., Christodoulou, Shlosman & Tohline 1995 and refs. 
therein), with the exception of Athanassoula (2002). It is 
known that rigid halos are less hospitable to bar growth. 
All equal, bar instability in this case requires larger disk- 
to-halo mass ratio than in the live halo systems. This be- 
havior of stellar disks embedded in non-responsive halo 
potentials can be readily understood in terms of absence 
of (resonance) halo orbits capable of absorbing the disk 
angular momentum. As expected, it is accompanied by 
a nearly constant bar pattern speed, as the angular mo- 
mentum loss from the bar is minimized to the outer disk 
only 

Fig. 121 exhibits some properties of stellar bars in such 
rigid axisymmetric (RS) halos — with a large flat-core 
i?H = 10.0kpc (model RS 1), a smaller core i?H = 2.0 kpc 
(model RS4) and a cuspy (model RS5; i?H = 0.5 kpc) 
density profile (see Table^ . The (normalized) amplitude 
of the m = 2 mode, A2 (Fig. |31 top panels) provides only 
a partial description of the bar evolution. Because it 
is an integrated quantity over some specific radius and 
because the bar sizes will differ from model to model, 
it does not allow for a direct comparison between bars 
in different models. One can encounter a small size but 
strong bar whose amplitude will be diluted (say) by a 
large averaging region. To resolve this dilemma, we use 
both the bar amplitude, measured within a thin cylinder 
of radius of 5 kpc, and the bar length and its maximum 
ellipticity, e = 1— b/a (obtained from the isodensity fits, 
in the plane (Fig. |21 middle panels) , to characterize its 
strength. Both approaches can be tested observationally. 

We note that a reliable determination of a iV-body bar 
size from isophote fitting has its pitfalls (see Martinez- 
Valpuesta & Shlosman 2004). Before the vertical buck- 
ling, the bars have a flat distribution of ellipticity with 
radius, while at later times the ellipticity has a clear max- 
imum. We find empirically that a good approximation to 
the bar size would be the radius where the bar ellipticity 
drops by 0.1. This is verified using nonlinear orbit anal- 
ysis, finding the x-extent of the largest stable xi orbit in 
the bar (Martinez- Valpuesta, Shlosman & Heller 2006). 
On the other hand, the ellipticity of the surrounding disk 
shows a gradual decline with radius, even in mildly tri- 
axial halos used in this work. This different behavior of 

^ We are interested only in this property of spherical halos and 
refer to them as axisymmetric 



ellipticity in numerical bars and disks allows us to safely 
separate them. 

The rigid axisymmetric halo models RS 1 , 4 and 5 have 
been arranged along the same sequence of decreasing core 
size as in Fig. El These models have the same disk-to- 
halo mass ratio within the central 10 kpc, but correspond 
to increased central concentration in the halo. 

All three models develop a strong bar in the process of 
a 'normal' bar instability (Fig. 3). The bar amplitudes, 
show a strong peak after the initial bar growth, at 
T ~ 30 — 50, and a subsequent decline, partly associated 
with a transient m = 3 mode A^. The latter appears 
to be a purely numerical artifact resulting from mix- 
ing of analytical (halo) and live (disk) potentials. The 
more centrally-concentrated models supress the (planar!) 
m = 3 mode more efficiently. The sudden weakening and 
shortening of the bars has been related to the onset of 
chaos in strong bars (Martinez- Valpuesta & Shlosman 
2004), which can be the case in these models as well. 
Asymptotically, the more cuspy models show marginally 
stronger bars. While the RS4 bar exhibits a slight secu- 
lar growth in A2, the RS 5 bar shows a secular decline at 
the same time, again probably related to stronger chaos 
excited in cuspy models. 

The asymptotic bar size, ~ 6 — 7 kpc, is largest in RS 1 
— the large halo core model, and about 3 — 4 kpc in RS5. 
The maximal bar size (at r ~ 30 — 50) also clearly cor- 
relates with the halo core size, i?H- This trend is much 
more pronounced in models which have increasing cen- 
tral mass concentration (see section 3.2.1 and Fig. 5). The 
bar ellipticity, e, closely follows the size evolution of the 
bars in all models. The bar pattern speed, lip, shows 
some initial decline with time, more substantial in cuspy 
models, but then levels off — as expected, because of lack 
of angular momentum absorbing material in axisymmet- 
ric and to a certain degree also in mildly triaxial halos. 
The outer disk quickly saturates for this redistribution 
of angular momentum. As expected, more cuspy models 
have progressively faster flp. 

We find that the bars in rigid halos are not subject to 
the buckling instability (e.g., Pfcnnigcr & Friedli 1991), 
apparently because of difficulty to develop vertical asym- 
metric modes, m = 1. But vertical m = 3 and 5 show up 
progressively more in RS 4 and RS 5, at the level of ~ 8%. 
The vertical resonant heating in the bar (e.g., Friedli & 
Pfcnnigcr 1990) is not obviously observed, possibly due 
to the rigid potential of the halo, and especially due to 
the absence of a discreteness noise in the analytical po- 
tential. 

The effect of the bar on the disk is to shorten the 
original radial scalelength in the inner (bar region) disk 
and in the outer disk. In this sense, the combination of 
initial conditions and disk evolution maintain a double- 
exponential disk, td, characterized by the inner and 
outer radial scalelengths. The inner disk (the bar region) 
in RS 1 develops td ~ 2.2 kpc, when measured along the 
bar major axis. This value is progressively smaller for 
more cuspy models, e.g., ~ 1.6 kpc in RS 5. In the outer 
disk, td decreases from its initial value 5.7 kpc to about 
5 kpc with time and remains constant when disk/bar 
reaches quasistationary stage after r ~ 80 (in RS 1). The 
Toomre's Q increases overall from 1.2 to 2.1, measured 
at two radial initial scalelengths after the bar weakening 
in RSI and more so in RS 4 and RS 5. 
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Fig. 3. — Bar strength (top panels; m = 2 amplitude A2: solid line, m = 3 amplitude A3: dotted line), bar length and (max.) ellipticity 
(middle panels; solid and dotted line, respectively) and pattern speed Hp (bottom panels) as a function of time for models RSI, 4 and 5 
(from left to right) wit rigid axisymmctric halos. 
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Fig. 4. — Same as Fig.|21 but for models with rigid triaxial halos, i.e. RT 1, 4 and 5 - constant mass ratio within 10 kpc. 



The vertical disk scaleheight, zd, shows progressively 
less change with increased cuspiness. While it increases 
from 0.5 kpc to ~ 0.7 kpc in RSI, as a result of disk 
heating, towards the end of the simulations, it increases 
only 0.6 kpc in RS 4 and stays unchanged in RS 5. 

3.1.2. Rigid triaxial halos 



In this section we describe the evolution of the models 
with rigid triaxial halos and a constant A/d /-^h within 
lOkpc (TablcUJ, i.e. models RT 1, RT4 and RT5. Some 
of the bar properties in these models arc shown in Fig.^ 
The angular momentum redistribution is discussed in 
section 4. 

The bar initial growth and size, its maximal strength 
(i.e., of A2), and even its peak ellipticity in rigid triax- 
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ial halos arc remarkably similar to those in axisymmet- 
ric ones. The only difference is that the growth rate is 
faster and so the peak strength is achieved earlier. The 
subsequent evolution, however, is very different — the 
bar dissolves over the characteristic time of 30 — 100, 
depending on the model, leaving only a weak oval distor- 
tion behind (Fig. 4) . The general trend of bar evolution 
in this sequence of models clearly links the halo cuspiness 
to the dissolution process, which is amplified by the halo 
triaxiality — a trend predicted by ES02 based on the or- 
bital stability analysis and the development of connected 
chaotic domains in these models 

The bar size plot (mid panels), which indicates a bar of 
6 kpc in RT 1 at later times, is misleading here and seems 
to represent only the disk response to the triaxial poten- 
tial and the ellipticity is that of the disk. As long as the 
bars exist in these models, their pattern speeds mirror 
those in axisymmetric halos of the same mass concentra- 
tion. 

Despite the destruction of bars in RT 1 to RT 5 models, 
the disk tends to acquire a double-valued radial expo- 
nential scalelength: the inner one of about 2 kpc and the 
outer one of ^ 5—6 kpc. The Toomre's Q and the vertical 
scaleheight behave as in the axisymmetric sequence. 

3.2. Comparing Bars in Rigid and Live Axisymmetric 
Halos with Decreasing Disk-to-Halo Mass Ratios 

3.2.1. Rigid axisymmetric halos 

To provide a more direct comparison to the results of 
ES02, we present in this section a set of simulations in 
which we closely follow their model parameters. The ba- 
sic difference with the previous set of models is that halo 
mass within the central 10 kpc is increased for progres- 
sively smaller cores (Table 1). Thus, when moving from 
flat core to cuspy models, we move from a maximal disk 
model (RS 1) to a halo-dominated one (RS3) as evident 
from Fig. 2. Fig. shows the evolution of stellar bars in 
such rigid axisymmetric halos. This corresponds to the 
increased central concentration, in excess of that shown 
by RS l^RS 4^RS 5 sequence. 

Fig. 5 exhibits the disk evolution inside such halos. 
The strong bar only develops in the model RSI. The RS2 
supresscs the bar growth for the first r 60, developing a 
short, ~ 2 — 3 kpc bar. Thereafter, A2 drops and the bar 
further weakens to a rather oval distortion which persists 
over the period of simulations. The RS3 supresses the 
bar growth rate even stronger, and A2 shows a steady 
secular growth not leveling off even after r ^ 300, while 
A3 is completely supressed in a more cuspy models. 

Q increases to ~ 2 while the vertical scalelength shows 
no change whatsoever. The bar pattern speed correlates 
with the halo's mass concentration, i.e. increasing with 
smaller core sizes. Pattern speed in RS 3 is linearly de- 
creasing with time reflecting the secular increase in A2 . 

To summarize the sequence of rigid axisymmetric halos 
RS 1— i'RS 3: the stellar bar is basically supressed in more 
cuspy models in accordance with ES02 analysis. 

3.2.2. Live axisymmetric halos 

Models LS 1 to LS 3 present the case of live axisym- 
metric and increasingly cuspy halos. After introducing 
the disk and letting the halo to relax, we have fitted the 
halo density profile at the end of this relaxation process 



using the parameters of logarithmic potential, Vh, ^h, 
f3 and 7 (see Table QJ. The new density profiles for the 
halo can be fit reasonably well and show slight increase 
in V}i and unchanged p. Owing to the flattened disk 
potential 7 slightly decreases. For model LS 1, the core 
radius remained at ~ 10 kpc, for model LS 2 it decreased 
to ~ 1.8 kpc, and for model LS3 — to some 0.3 kpc. 
The main difference between these models and the pre- 
viously discussed ones with rigid halos, however, is an 
active redistribution of angular momentum between the 
bar forming region in the disk and the live halo (see sec- 
tion 4). For all live axisymmetric halos during disk evo- 
lution the halo density proflle is stable (e.g., as shown for 
LS3 in Fig. 7). 

The sequence LS1^LS3 has produced the strongest 
bars among all of our models (Fig. 6). The amplitude 
A2 has also shown a pattern of behavior: (1) the bar 
growth is fastest in LS 1 (disk dominated) and slowest in 
LS3 (halo dominated); (2) the maximal bar strength is 
increasing along the sequence; (3) all bars exhibit ver- 
tical buckling, (4) A2 decays most strongly in LSI af- 
ter the buckling, and least in LS 3; and flnally (5) the 
post-buckling bar in LS 1 continuous its secular growth, 
which quickly saturates in LS2 and starts to decay in 
LS3. Hence LSI (and to a lesser degree LS2) shows 
most resemblance to the secular evolution of a stellar bar 
in Martinez- Valpuesta, Shlosman & Heller (2006) where 
the bar weakens dramatically during its buckling and in- 
creases its strength thereafter. 

For model LSI, the bar size decreases from roughly 
10 kpc just before the buckling, to ~ 7 kpc at r ^ 60 
and then increases again to about 9 kpc after the buck- 
ling. For model LS 3, it grows from 5 kpc to 8 kpc. This 
means that during the initial growth, up to the time of 
the buckling, the bar is conflncd to the regular region 
delineated in ES02. 

The bar sizes show clear correlation with the halo core 
size, i?H, as in the rigid sequence RSI RS5 described 
in section 3.1. One can clearly observe here the 'delay' 
in the bar growth in halo-dominated models. The bar 
ellipticity shows a visible decay after the buckling in LSI 
and stabilizes afterwards, while in LS 2 and LS3 e does 
not decay. 

As expected, the bars slow down substantially over the 
simulation time and the initial slowdown is well corre- 
lated with i?H (Fig. El bottom panels). This angular 
momentum transfer from the disk is deposited mostly in 
the internal circulation in the halo as its figure docs not 
acquire rotation (section 4). 

We have mentioned above that all LS models show 
buckling which is less pronounced for halo-dominated 
models. This can be noted from the amplitudes of the 
vertical modes — LS 1 shows clear increase of to = 1, and 
LS 2 and 3 only in m = 3 and 5. 

The disk evolution in LS models again leads to a 
double-exponential scalelength: the inner bar-dominated 
part has td which slightly decreases from LSI to LS3, 
from ~ 2 kpc to ^ 1.6 kpc, while the outer disk scale- 
length actually increases along the sequence, from about 
5 kpc to just below 6 kpc. 

The disk scaleheight in LSI shows an abrupt increase 
at T^50 from its initial value 0.5 to about 0.9 kpc within 
few dynamical times, during the vertical buckling of the 
bar. This effect is clearly related to the formation of a 
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Fig. 5. — Same as Fig.U] but for models with rigid axisymmetric halos, i.e. RS 1, 2 and 3 - from disk to halo dominated. 
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Fig. 6. — Same as Fig|3] but for models LS 1-3 with live axisymmetric halos. 



pseudo-bulge, i.e., bulge of a boxy/peanut shape. LS2 
experiences a milder increase in zd to ~ 0.65 kpc and in 
LS 3 it remains largely unchanged. 

All LS models show the formation of pseudo-bulges as- 
sociated with the vertical buckling of the bars. However 
the shapes of these bulges differ and this difference per- 
sists for the time of the simulations. While in model 
LS 1 and LS 2 the isodensity contours show a more boxy 
shape, the cuspy model, LS 3 develops a peanut-shaped 



bulge (see also Athanassoula & Misiriotis 2002). 

3.3. Comparing Bars in Rigid and Live Triaxial Halos 

3.3. L Rigid triaxial halos 

The RT1^RT3 sequence describes the disk evolu- 
tion within rigid triaxial halos with decreasing core size, 
closely following the initial conditions in ES02 (Fig. 8). 
The only difference between the previously discussed 
RT1^RT5 sequence and RT 1-3 is that the latter one 
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Fig. 7. — Density profile along the halo main (longest) axis for 
time t = (thin line) and t = 300 (thick line) for model LS 3 
(rn = 0.05). 



is progressively halo-dominated. This difference is re- 
sponsible for the damping of bar instability in a disk 
embedded in RT 2 and RT 3 halos. Even RT 1 shows a 
substantial decay of the initially strong bar to a largely 
an oval distortion, unlike the bar in RS 1 model. Hence, 
asymptotically RT 2 and RT 3 behave as RT 4 and RT 5. 
However, early in the disk evolution, the former models 
show strong bars (see section 3.2.1). 

The size and cUipticity of the central (bar) oval distor- 
tion in RTl to RT3 correlate well with Ru. In RTl, 
the maximal bar length is ~ 9 kpc at t ~ 35 and after 
T ~ 50 — 60, the bar deteriorates into a rather triaxial 
(in the equatorial plane) configuration within the central 
r ~ 4 — 5 kpc. Note the A2 amplitude becomes finite 
already in the first moments into the simulation because 
of the disk response to the halo triaxiality. 

The first few disk rotations, especially in RT2 and 
RT3, show a grand design spiral structure in the outer 
disk which remains tightly bound and quickly decays. 
Even more than in RS3, the disk stays largely axisym- 
mctric within the central few kpc after r ~ 30. Its sur- 
face density profile shows (almost) no evolution. The 
outer disk radial scalclcngth decreases to about 4.5- 
5 kpc, while the disk scalchcight stays unchanged. 

3.3.2. Live triaxial halos 

As in the case of the live axisymmetric halos, the tri- 
axial halos are modified by the introduction of the frozen 
disk potential. At the moment of disk 'release,' the 
halo core size in model LT 1 has slightly increased to 
i?H ~ 10.5 kpc, in LT2 it has increased to ^ 2.4 kpc, 
and in LT3 stayed at 0.5 kpc (see Tabic 1). We found it 
increasingly difficult to maintain the (inner) halo triaxi- 
ality in cuspy models after introduction of a frozen disk. 
This has affected the LT 3 model which possesses a more 
structurally unstable halo. 

The evolution of A2 amplitude in the LTl— >LT3 se- 
quence is remarkably different from that in the rigid se- 
quence RT 1 to RT 3. The initial disk response to the 
halo potential has induced strong bars in all these mod- 
els. The fate of the bar, however, differs depending on 
the model. The bars decay nearly linearly in LT 1 and 
LT2 on a timescale of t ~ 100, and enter a steady state 
in LT 3. Overall, the LT 3 model behaves differently and 
its evolution resembes more the LS 3 model, although the 
bar is not so strong. The buckling is clearly visible in all 
LT models. However, the resulting boxy shapes of bulges 
are quickly erased in LT 1 and LT 2, while peanut shape 



bulge in LT 3 persists till the end of the simulation. The 
planar A3 amplitude is negligible in all live halo models. 

The radial density profile in the triaxial halo is stable 
during live disk evolution in all these models, even for the 
cuspy LT3. However the halo triaxiality T is reduced 
sharply in LT3, due to decrease in en- After t ~ 70, 
it has largely lost its prolateness, (3^1 (Fig. 10). This 
halo axisymmetrization is only partially a result of initial 
conditions — LT 3 has had somewhat larger /? at t = 0, 
as noted above. Rather, as we discuss in section 5, the 
halo is much more structurally unstable in cuspy models. 
In order to test this latter assumption, we have re-run 
the LT2 model with a frozen disk. Fig. 11 shows the 
resulting evolution of potential axial ratios for the halo. 
Both halo flatness, /, and its prolateness, en, show no 
evolution over 300 dynamical times, ~ 14 Gyrs. This 
confirms that it is the developing bar perturbation which 
decreases the halo prolateness without much effect on its 
flatness. 

The bar size and ellipticity LT 1 and LT 2 (Fig. 9) rep- 
resent a combined disk and bar evolution. The outer 
disk develops elaborate system of long-lived spiral arms, 
much more pronounced than in live axisymmetric mod- 
els. Because of the difficulty to disentangle between the 
bar, spiral arms and the oval disk, the estimated bar size 
in Fig. 9 is erroneous. The inner disk radial scalelength 
shows a general decrease to about 1.7 kpc in LT 1 and to 
~ 1.8 kpc in LT2 and LT3. The outer disk has td ap- 
proaching ~ 4.8 kpc in LT 1 and about 5 — 5.5 kpc in LT2 
and LT3. The inner disk scaleheight increases abruptly 
to ~ 0.9 kpc after r ~ 35 in LT 1 and to ^ 0.8 kpc in 
LT2 and ~ 0.75 kpc LT3 models. 

The live triaxial models are especially efficient in trig- 
gering the spiral structure in the barred disk. The arms 
penetrate deeper towards the center in LT 2 model with 
smaller core radius. The strength of the spirals depends 
on the position angle of the bar with respect to the 
longest halo axis — when the bar is normal to this axis, 
the spirals are prominent. This type of spiral regenera- 
tion happens as long as the bar exists. 

3.3.3. Models with continuous support for halo 
triaxiality 

Models LT 1 to LT 3 show that the disk (and bar) re- 
sponse to a live triaxial (halo) potential can act as to 
decrease its prolateness, e.g., equatorial ellipticity. This 
effect is most pronounced in model LT 3, where it has 
been difficult to induce the halo's triaxiality, especially 
in its central region. Compared to the models with larger 
cores, the disk in LT 3 evolves more similarly to that in 
axisymmetric halos. Because of the strong feedback of 
stellar bars onto the halo triaxiality, we have run a model 
for LT 3 with a continuous support for triaxiality — the 
LT3MT model. The simple justification for this is that 
mergers will induce triaxiality in the extended DM ha- 
los in the first place. We, therefore, mimic this process 
by maintaining or continuously regenerating a mild halo 
triaxiality, and at the same time allowing the live halo to 
interact with the live disk. 

Fig. 12 (right) exhibits such an evolution of LT 3 MT 
and provides A2 from LT 3 for a comparison. The adi- 
abatic heating/cooling procedure was gradually turned 
on and has its full strength from r = 70 to 225, when 
it is gradually turned off. Initially, the LT 3 MT evolved 
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Fig. 8. — Same as Fig.U] but for models with rigid triaxial halos. 
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Fig. 9. — Same as Fig. 13 but for models with live triaxial halos. 



identically to LT 3. Because the halo triaxiality has been 
maintained at the steady level (Table 1), it has a pro- 
found effect on the disk. The bar completely dissolves 
and starts to grow afterwards when the triaxiality is 
gradually washed out. 

3.3.4. Light disks in triaxial halos 

The interaction between the halo and the evolving disk 
tends to gradually wash out the initial triaxiality in the 



halo. The efficiency of this process depends on the disk- 
to-halo mass ratio. The bar is then able to grow within 
the region in which the halo becomes more axisymmet- 
ric. To show explicitly the dependency of halo and disk 
evolution on their (inner) mass ratios, we have run live 
triaxial models with lower disk-to-halo mass ratio within 
central 10 kpc for the cuspy LT3 model. Here we show 
only LT3HM — the model with a half disk mass. The 
bar instability in this model is substantially suppressed 



12 



0.8 
0.6 
0.4 
0.2 





LT 1 



Berentzcn, Shlosman and Jogee 

LT 2 




LT 3 



Fig. 10. — Halo (isopotcntial) axes ratios /3 (full lines) and 7 (dotted lines) as a function of radius for models with live triaxial halos. 
Thin and thick lines indicate the profiles at the beginning and the end {t= 300) of the runs, respectively. 
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Fig. 11. — Testing the halo evolution in the LT2 model: re-run 
of the model with a frozen instead of the live disk, for comparison. 
Axes-ratio /3 (full line) and 7 (dotted line). The thin and thick lines 
represent the times t = and t = 300, respectively. The absence 
of evolution in f3 and 7 in this case confirms that the evolution 
observed in Fig. 10 results from interaction between the disk and 
the halo. 



compared to LT 3 model. LT 3 HM develops a weak oval 
distortion in the center (Fig. 12, left) whose shape de- 
pends strongly on its orientation with respect to the halo 
major axis — its ellipticity is smaller when it is normal 
to the halo, and larger, when it is parallel. Otherwise, 
the disk is stable and no bar instability develops. 

4. DISK-HALO ANGULAR MOMENTUM EXCHANGE 

The general direction of the angular momentum flow in 
all models is expected to be from the inner bar-forming 
disk towards the outer disk and the halo. This reflects 
behavior of all models that break their axial symmetry 
for whatever reason, that of initial conditions (triaxial 
halos) or as a result of a spontaneous bar instability. 
However, when this redistribution is mediated by the bar, 
there is also an additional flow from the inner disk (inside 
the bar corotation) to the outer disk. 

In the models with rigid axisymmetric halos the ex- 
change of angular momentum is limited to that between 
the inner and outer disks, across the corotation and me- 
diated by the bar. Models RS 1 and RS 5 (Fig. 13) show 
that the total loss of the angular momentum from the 
disk is ALz = and that the outer disk disk gains the 
angular momentum from inner part. This gain is larger 
for RS 5 as it develops stronger bar. Rigid triaxial mod- 
els drain some of the angular momentum from the disk 
(e.g., RT5 in Fig. 13) but this does not dominate the 
L-exchange because of the mild triaxiality here. 

Fig. 14 shows the effect of introduction of a live halo, 
axisymmetric and triaxial. The total loss of the L by 
the disk has been dramatically amplified compared to 
the rigid triaxial halos. The outer disk still gains the 



same amount but most of L goes to the halo. The strong 
and long-lived bar developing in the LS 1 model is in- 
strumental in this transfer. The main difference with 
the LTl model is that the disk immediately acquires a 
strongly oval shape which is gradually lost and so overall 
the transfer amount of angular momentum is less than in 
LSI. However, during the first 20-30 dynamical times, 
the L-transfer is stronger in LT 1. Subsequently, the bar 
which has developed in LT 1 dissolves, which levels off 
the rate of angular momentum trasfer. 

Furthermore, in order to test the effect of the resonance 
coupling between the disk/bar and the halo, we have per- 
formed the spectral analysis of the orbital motion (Bin- 
ney & Spergel 1982) in conjunction with our nonlinear or- 
bital finder algorithm (Heller & Shlosman 1996). We find 
that a fraction of disk particles in live axisymmetric halos 
is locked by the lower resonances, especially by the inner 
Lindblad resonance (ILR). A fraction of the halo particles 
is locked by the corotation resonance — hence ILR-CR 
resonance coupling. Following Athanassoula (2002), we 
have frozen the model potentials at times t = 100 and 
150, and compared the angular momentum of individual 
particles at these times. Most of the angular momentum 
flow between the disk and the halo was mediated by this 
ILR-CR coupling. Particles trapped at t = 100 largely 
remained trapped at 150. This result shows that the 
adverse effect of the grainy iV-body potential on the evo- 
lution in our models is limited, and the latter is driven 
by a chaotic dynamics instead. 

5. DISCUSSION AND CONCLUSIONS 

We have analyzed the formation and evolution of stel- 
lar bars in galactic disks embedded in rigid (i.e., analyti- 
cal) and live axisymmetric and mildly-triaxial DM halos 
of a varying cuspiness. The nearly axisymmetric poten- 
tial employed in our simulations allows us to focus on 
the development of bar instability per se and does not 
violate the overall dynamics and survival of the disks. 
Using tailored numerical simulations we aimed at testing 
and extending the predictions of El-Zant & Shlosman 
(2002, ES02) which employed the method of Liapunov 
exponents to address the fate of bars in analytical tri- 
axial halos. Our simulations of live halos, unlike ES02, 
include the feedback between the disk, bar, and halo. We 
fully support the main conclusions of ES02 and provide 
some additional insight into the bar and disk evolution. 

We start with summarizing our results and follow up 
with discussion — first based on the bar strength. In 
all cases we can separate the initial bar growth to its 
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Fig. 12. — The bar amplitude A2 for models LT3HM (0.5 kpc core) with a half mass disk (left panel) and LT3MT continuous support 
of triaxiality (right panel). The dotted line shows A2 of LT3 for comparison. 
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Fig. 13. — Models RSI, RS 5 and RT5 showing the angular momentum exchange in the live disk along the sequence with increasing 
cuspiness and with a rigid axisymmetric (RS 1, RS 5) and triaxial (RT 5) halos. The change in the angular momentum in the disk, AL^ , is 
normalized by the total angular momentum in the disk, Lz, Oi at r = when it is released. The solid line represents integrated ALz / L^^a 
in the disk between r = — 10 kpc, the dotted line — for r > 10 kpc and the dashed line — the change in the total angular momentum in 
the disk, which is absorbed by the halo. 
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Fig. 14. — Models RS 1, LS 1 and LT 1 showing the angular momentum exchange in a rigid vs live axisymmetric vs live triaxial halos. 
As in Fig. 13. 



full strength, over r ~ 30 — 50 ~ 1.5 — 2 Gyr, from the 
subsequent bar evolution, dynamical or secular, over the 
Hubble time. The first stage is very similar between all 
rigid axisymmetric and triaxial halo models. We find, 
that while some of the rigid halos suppress the bar for- 
mation, depending on the disk-to-halo mass ratio D/H 
within the disk radius, all live halos develop impressive 
bars with ellipticities in the range of ^ 0.4 — 0.7, de- 
pending on their evolutionary stage. In axisymmetric 
live cases (only!) these bars appear stronger for smaller 
D/H, in line with Athanassoula & Misiriotis (2002). In 
live halos, the bars develop faster in triaxial compared to 
axisymmetric halos. (This statement is marginally true 
also for rigid halos.) 

All bars weaken at the end of the first stage. This 
process is accompanied by bar's vertical buckling in the 
live halos and formation of exponential (pseudo-) bulges 
of a peanut or boxy shape. No buckling is observed in 



rigid halos. 

The subsequent evolution of the bar differs substan- 
tially between axisymmetric and triaxial halos. In the 
former case, the bars appear dynamically stable and 
show a limited secular evolution by a factor of 2 in 
strength, in growth or decay after the initial vertical 
buckling. This result is in agreement with Athanassoula 
(2002). In the latter case, the bars dissolve almost com- 
pletely on a timescale of r ~ 30 — 100 ~ 1.5 — 5 Gyr, 
sometimes leaving behind a weak oval distortion in the 
central regions. This is true for all models with the ex- 
ception of LT 3 — a live cuspy triaxial halo, where the 
triaxiality (i.e. prolateness) is erased early due to the 
feedback from the bar and its ensuing evolution closely 
resembles that of the axisymmetric halos. 

To verify that the bar endurance is related to the era- 
sure of the halo prolateness, we have run a test model 
where (live) halo triaxiality in LT 3 has been maintained 
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at a constant level. The bar dissolved in this case as ex- 
pected. On the other hand, to show that the evolution 
of halo prolateness is due to destabilizing action of the 
bar, we run models with a frozen axisymmetric disk in 
live triaxial halos — these models show basically no evo- 
lution for Efj- We have also run a number of test models 
with the LT3 halo but smaller D/H mass ratio. Mod- 
els with D/H^ 0.5, formed only a weak oval distortion. 
This is unlike the trend in the axisymmetric live halos, 
where small D/H resulted in stronger bars. We discuss 
this below. 

Second, the bar sizes appear to correlate with the halo 
core sizes, with the exception of live triaxial models 
where the bar size is more difficult to establish due to 
the strong response in the disk and a quick bar disso- 
lution. The bar ellipticities provide additional source of 
information on the bar strength and even more about the 
state of the disk when the bar amplitudes measured by 
A2 decay below 0.1. 

Third, the pattern speeds of bars anticorrelate with 
halo core sizes — more centrally concentrated halos pro- 
duce faster tumbling bars. However, the deceleration 
rate of the bar anticorrelates with the halo core sizes 
as well — thus, after the Hubble time, the end pattern 
speed of the bars is similar. As expected, bars in rigid 
halos quickly reach constant pattern speeds, while bars 
in live halos continue to slow down asymptotically. In- 
terestingly, bars in live axisymmetric halos can become 
stronger or weaker while slowing down. Bars in live triax- 
ial halos always weaken and dissolve while slowing down. 

Fourth, in rigid halos, the angular momentum of the 
inner disk is fed into the outer disk, across the bar coro- 
tation. Some angular momentum is lost to the halo in 
rigid triaxial models. In live models, the halo appears to 
gain angular momentum, more so fractionally in triaxial 
models. We have verified that in axisymmetric models 
the angular momentum is mediated by lower resonances, 
especially by the inner Lindblad resonance in the disk 
and by the corotation resonance in the halo. 

Finally, in triaxial models, as long as a bar is present, 
the disk evolution is characterized by development of a 
strong spiral structure outside the bar region. It is es- 
pecially prominent in live triaxial halos. The outer disk 
acquires exponential surface density distribution, while 
the inner (bar-unstable) part can be modeled roughly by 
a shorter exponential scalelength when measured along 
the bar major axis. 

Our emerging understanding of the evolution of bars 
embedded in live, mildly triaxial or axisymmetric halos 
is based on two physical processes which determine the 
fate of the system. These are the angular momentum 
redistribution in the disk-halo system and the develop- 
ment of chaos. We can make a general statement that 
angular momentum flows from the bar unstable region 
in the disk to the outer disk and the halo. This pro- 
cess is accompanied by the initial growth of the bar both 
in strength and in size. In axisymmetric halos, it ulti- 
mately saturates and the bar enters a steady-state phase 
characterized by secular changes (Athanassoula 2003). 
Instead, in a mildly triaxial halo which hosts a barred 
disk, the dominant process is the increase in the frac- 
tion of chaotic trajectories — it affects both the bar and 
halo structures. In other words, no steady-state devel- 



ops, and either the bar is dissolved or halo prolateness is 
washed out. This effect on the bar embedded in a rigid 
triaxial halo has been calculated in ES02, and here we 
have presented fully self-consistent numerical simulations 
accounting for the bar and live halo evolution. 

After the submission of this work we have learned 
about the work by Curir et al. (2005) where stellar disks 
have been evolved within the cosmological DM halos, al- 
beit at much lower resolution than in our work. In these 
simulations, the bars form in responce to the original 
torquing by the halo, but, contrary to our main results, 
survive for the rest of the simulation. We strongly sus- 
pect that this different behavior of bars results from the 
loss of prolateness by DM halos in their models, simi- 
larly to our model LT3. Furthermore, it is the process 
of introduction of disk into the halo which is also respon- 
sible for washing out the halo prolateness (Section 2.1), 
which required us to apply the second heating-cooling 
procedure. Unfortunately, Curir et al. do not provide 
any analysis of the shape of their halos after bringing up 
the disks and during their subsequent evolution. Their 
Table 1 refers to purely DM halos without stellar disks 
only. 

Thus we expect that the combination of halo triaxiality 
and its cuspiness in conjunction with a fast rotating^ bar 
will lead to the development of chaos and to the disso- 
lution of the least massive or structurally stable triaxial 
feature in the system, i.e., the bar or halo triaxiality. 
In section 1 (see also ES02 and refs. therein), we have 
argued that cuspy density distribution cause solutions 
for the Poisson equation to be far from quadratic and 
this will produce the coupling between different degrees 
of freedom in equations of motion which become highly 
nonlinear — a prime recipe to develop a chaotic behavior 
when such a system is perturbed. 

To understand the results of these numerical simula- 
tions within the context of developing chaos in a triax- 
ial system with a fast tumbling stellar bar, we take the 
next step and discuss the evolution of the bar strength in 
terms of the stability of trajectories used by ES02. Ul- 
timately bars and other morphological features, such as 
prolate halos, dissolve when the density figure stops to 
support the figure of the background potential, in other 
words when self-consistency of density-potential pair is 
violated. As mentioned in section 1, a stable 3-D figure 
must be built from trajectories which at least approxi- 
mately conserve the invariants of motion. Chaos appears 
when the number of invariants of motion becomes smaller 
that the dimensionality of the system. 

Two issues appear most relevant in this respect: the 
overall asymmetry of the background potential and the 
degree of mass concentration. Stability of trajectories 
for potentials and density distributions used in this work 
have been analyzed in ES02 by calculating the Lia- 
punov exponents on a high-resolution 3-D cylindrical 
grid. These exponents provide a measure of timescales 
which are associated with dynamical instabilities and are 
especially suitable to study development of chaotic mo- 
tions in time-dependent potentials (see detailed descrip- 
tion of this method in sections 1, 2 and the Appendix of 
ES02). 

The mapping of halos which arc under consideration 

^ As defined in the footnote of section 1 
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here using the Liapunov exponents has shown that tri- 
axial halos are regular in the absence of a bar — the top 
panel of Fig. 6 in ES02 indicates that trajectories within 
this potential are stable well over the Hubble time. This 
is supported by live models of triaxial halos presented 
here. All the models appear stable unless bar-unstable 
disks are added, even the cuspy halo in LT3, as expected. 
It is the addition of a (triaxial) bar fast tumbling with 
respect to the halo, which destabilizes the trajectories 
and triggers the chaos. The middle panels in Fig. 6 of 
ES02 show this effect in a cuspy halo — most of the 
configuration space becomes chaotic with the character- 
istic timescale ^ 1 Gyr. Moreover, the chaotic regions 
are fully connected, thus providing a strong indication 
that survival of any non-axisymmctric structure is highly 
questionable. The less cuspy barred models of ES02 show 
progressively less chaotic regions, especially in the flat 
cores, and the existing chaos is less interconnected. A 
self-consistent bar-disk model of Pfenniger (1984), with- 
out a surrounding halo, shows chaotic regions confined to 
the bar corotation radius (lower panels. Fig. 6 of ES02), 
in line with the above explanation. We also note that the 
remarkable correspondence between the distributions of 
the available configuration space volume for the trajecto- 
ries in ES02 and the maximal Liapunov exponents points 
out that the nearby stability islands are not efficient in 
trapping the chaotic trajectories. Therefore, one should 
expect that the predicted bar dissolution in ES02 models 
will indeed be observed in numerical simulations. 

The sequences, RTl^RTS and RT1^RT5 discussed 
in section 3 represent a gradual increase in the halo cuspi- 
ness, i.e., increased central mass concentration and pro- 
gressively smaller size of the halo core. Figs. 3 and 6 
of ES02 demonstrate that the regular (stable) region 
shrinks from ~ 5 — 6 kpc in RT 1 to about 0.5 kpc in 
RT3, along with the increased halo cuspincss, not leav- 
ing much space for the bar to survive. When the bar 
feedback onto the halo is fully accounted for, i.e., in live 
triaxial halos — they also lose part of their asymmetry 
(prolateness) en, but most of it remains. This explains 
why bars in LT 1 and LT 2 models behave qualitatively 
similarly to those in rigid RT models — all of them ex- 
hibit the bar dissolution. 

The main difference between these models and the 
cuspy LT 3 is that cuspy halos are more structurally un- 
stable, as we discussed above, and their prolateness is 
washed out early enough to allow for the bar to develop 
in a nearly axisymmetric environment. We confirm that 
either the onset of chaos dissolves the bar (e.g., LT 1 and 
LT2) or it destroys the halo prolateness (e.g., LT3). We 
strongly suspect that this latter process of is also at work 
in Curir et al. (2005) simulations of cosmological halos, 
allowing for long-lived bars. Therefore, to some extent 
long-lived bars are incompatible with the high (equato- 
rial) asymmetry in the DM halos. Even the mild asym- 
metry observed in some halos of nearby galaxies is suffi- 
cient to shorten the bar lifetime to ^ 5 Gyrs. 

This provides an interesting constraint on halo shapes, 
when taken in tandem with recent observational results. 
Using _ff5r-based morphologies and accurate rcdshifts 
from the Galaxy Evolution from Morphology and SEDs 
(GEMS; Rix et al. 2004) survey, Jogee et al. (2004, 
hereafter J04; 2005) recently showed that the optical 
fraction and distribution of structural properties for bars 



with moderate-to- high strength (e 0.35) remain similar 
from the present-day, out to look-back times of 2 — 8 Gyr 
(z ^ 0.2 — 1.0). J04 argue that these findings imply 
that on average bars have a long lifetime, well in excess 
of 2 Gyr. A constant optical bar fraction out to z 1 
is also reported independently from a smaller survey by 
Elmegreen et al. (2004). The simulations in this pa- 
per, when combined with these empirical results on a 
relatively constant fraction of bars out to z ~ 1 and an 
inferred long bar lifetime, put a lower limit on the halo 
equatorial axial ratio (3 = b/a ol 0.9 (in potential axes) 
at the redshift range of z ~ — 1. This corresponds ap- 
proximately to {b/a)p 0.75 — 0.8 axial ratio in density 
distribution. 

The disappearance of triaxial halos with en ^ 0.1 in 
disk galaxies at z ^ 1 seems as a corollary to our nu- 
merical simulations, when supplemented with the above 
observational results. As discussed in section 1, while 
cosmological simulations of dissipationless CDM galactic 
halos invariably produce triaxial halos, the triaxiality of 
the halo can be subsequently diluted by the addition of 
baryonic components. This is likely to happen, for in- 
stance, during the early formation and development of a 
galactic disk. Furthermore, this trend will be supported 
by a general decrease in the frequency of galaxy interac- 
tions and mergers below z 2. 

We note the following caveat: the halo prolateness in 
our modeled halos is (partially or fully) washed out by 
the action of the bar which leads to a dramatic increase 
in the fraction of chaotic orbits — these cannot support 
triaxial figures. It is the inner bar-forming part of the 
system where chaotic orbits dominate in asymmetric ha- 
los. In principle, we can envision the situation when 
the halo figure tumbles and trajectories which originate 
in its outer part cannot penetrate deep enough toward 
the central regions where they are typically destabilized. 
Such a halo will be more stable in preserving its prolate- 
ness, at least in its outer part. It is unclear how fast the 
halo must be tumbling for this effect to take place. Such 
models are beyond the scope of this work. 

Finally, we find it important to address the question 
to what extent the evolution of stellar bars in our A^- 
body potentials of mildly triaxial halos is not a numer- 
ical artifact, i.e., not a process driven by a numerical 
diffusion. In other words, can a Poisson noise, associ- 
ated with the discreteness of the A^-body potential, be 
primarily responsible for the observed model bar disso- 
lution, as opposed to the chaotic behavior triggered by 
the competing forces from the triaxial halo and a bar 
in ES02 for smooth analytical potentials? One can ar- 
gue that that the graininess of the A^-body potential is 
identical in the axisymmetric and triaxial halo models. 
Because there is no indication that this potential has 
caused any numerical 'relaxation' effects in the former 
models, the latter ones will not be dominated by these 
effects as well. In principle, the consequences of the dis- 
creteness noise in the 'mixed' (triaxial halo/bar) models, 
i.e., models which are intrinsically more chaotic, can be 
more complex than in the axisymmetric models. Never- 
theless, albeit indirectly, this supports our point of view 
that it is the chaotic dynamics which drives the evolution 
of stellar bar in triaxial halos modeled here. 

To provide the most direct answer to this question is 
to analyze the numerical models presented here in terms 
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of Lyapunov exponents of ES02, but this is outside the 
scope of our work and will be addressed elsewhere. How- 
ever, even if this approach will be taken, it is not without 
caveats as we deal here with a transient and not a clas- 
sical chaos which is defined in the asymptotic limit and 
corresponds to an exponential divergence of a trajectory 
in a /ixed potential (e.g., Lichtenberg & Lieberman 1995). 

In summary, we find that even a mild halo triaxiality 
of ~ 0.9 in potential axial ratio, dissolves the stellar bar 
on a timescale of ^ 5 Gyr, sometimes leaving behind a 
weak oval distortion. Especially in the live triaxial ha- 
los, as long as a bar is present, strong spiral structure 
develop outside the bar and its strength depends on the 
mutual orientation of bars and the halo major axis. In 
comparable live axisymmetric halos, there is limited sec- 
ular evolution, either growth or decay, of the embedded 
bars. In these cases, the halo core sizes correlate directly 
with the bar sizes and their central mass concentration 
— with the bar pattern speeds. Cuspy halos are more 
susceptible to washing out of their triaxiality due to the 
action of the bar, and the subsequent evolution is sim- 
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